home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 1 / Cream of the Crop 1.iso / EDUCATE / DIVECOMP.ARJ / DIVECOMP.SH / divecomp.c < prev    next >
Text File  |  1990-12-16  |  29KB  |  1,073 lines
  1. /* $Header: divecomp.c,v 2.5 90/12/12 15:30:02 dave Exp $
  2.    $Log:    divecomp.c,v $
  3.  * Revision 2.5  90/12/12  15:30:02  15:30:02  dave (Dave Waller)
  4.  * See file "rev_history" for information on new features.
  5.  * 
  6.  * Revision 2.4  90/11/27  15:30:03  15:30:03  dave (Dave Waller)
  7.  * new features:
  8.  * 
  9.  * - tissue loading graph now compresses when tissue loading exceeds
  10.  *   100%, to a scale of 0-200%. Rescales to 0-100% when all compartments
  11.  *   are <= 100%.
  12.  * 
  13.  * - Compartments can now be displayed as absolute pressure in (units) sea
  14.  *   water, instead of % loading. Default is loading graph, absolute
  15.  *   pressure can be selected at runtime with -P option. While the program
  16.  *   is running, user can toggle between the two modes by typing 't' (this
  17.  *   only works while calculations are taking place, not while waiting for
  18.  *   input). Scale along the top is calibrated in (units) sea water absolute.
  19.  * 
  20.  * - During calculation, depth can be adjusted up or down using the up and
  21.  *   down arrow keys on the keyboard. For example, if the program is calculating
  22.  *   a dive level of 60 feet for 50 minutes, the depth can be adjusted up or
  23.  *   down by pressing the arrow keys. Adjustments are made in increments
  24.  *   equal to the depth increment in the dive profile display.
  25.  * 
  26.  * - If compartments are overloaded, the program can "autodecompress" by
  27.  *   entering a depth of 'd' in interactive mode. The program will then
  28.  *   move up to the ceiling value and "ride" it until all the compartments
  29.  *   are <=100% loading.
  30.  * 
  31.  * - When displaying compartment loading, pressure values are in (units) SWA
  32.  *   instead of FSWA. (Of course, if the units are feet, then the values
  33.  *   *are* FSWA).
  34.  * 
  35.  * - Stats line during overloading now displays a more meaningful message
  36.  *   regarding ceiling and decompression time.
  37.  * 
  38.  * Revision 2.1  90/11/21  11:33:57  11:33:57  dave (Dave Waller)
  39.  * Changed SUN ifdefs to be macro SIMPLE instead... Makefile has been
  40.  * modified as well. Since Sun users can now compile with full SVID
  41.  * curses functionality, it seemed unnecessary to remove the reverse
  42.  * video bells and whistles for for them. However, this functionality
  43.  * is retained as a "simple" version of the program for those that do not
  44.  * have reverse video and underline capabilities on their terminals.
  45.  * 
  46.  * Revision 2.0  90/11/21  11:11:30  11:11:30  dave (Dave Waller)
  47.  * Bug fixes:
  48.  * 
  49.  * - incompatible typecasting in several places made the program fail
  50.  *   with certain model files.
  51.  * 
  52.  * - Logical error in the sample period determination algorithm.
  53.  * 
  54.  * Enhancements:
  55.  * 
  56.  * - added capability to display depth in arbitrary units. Default is
  57.  *   feet. Units are specified in either a profile file on the first line
  58.  *   or with the -u option. To use a different unit system, the user must
  59.  *   supply the unit name and a unit conversion factor that represents
  60.  *   units/ft. For example, to display in fathoms, the user would type
  61.  * 
  62.  *   divecomp -mEdge -ufathoms:0.1667
  63.  * 
  64.  *   There are 6 feet in a fathom, or 1/6=0.1667 fathoms per foot. For meters,
  65.  *   the specification would be -umeters:0.3048, or -um:0.3, etc.
  66.  * 
  67.  *   The unit specification is written into an output profile, so that the
  68.  *   correct units are used when using the profile. Units specified in a
  69.  *   profile file override command line specification.
  70.  * 
  71.  * - added logging capability. A log file can be specified with the -l
  72.  *   option, and the program will write the allowed nitrogen %age for
  73.  *   each compartment at each sample period into the log file. A header
  74.  *   containing information regarding the model, sample period, and depth
  75.  *   unit used is initially written into the file before the simulation
  76.  *   begins. To make a log of a computer run, type
  77.  * 
  78.  *   divecomp -mBuhlman -umeters:0.3 -llog.buhlman
  79.  * 
  80.  * 
  81.  * Revision 1.9  90/11/20  11:24:48  11:24:48  dave (Dave Waller)
  82.  * Bug fixes:
  83.  *      - logical errors in reading sample period from specified sources
  84.  *      - update_dive_profile() routine hada rounding error in the
  85.  *        code that compresses the display.
  86.  * 
  87.  * Enhancements:
  88.  * 
  89.  *      - new profile file format. Profiles are now stored as depth-duration
  90.  *        pairs, instead of discrete samples. This eliminates the need for
  91.  *        sample period specification in the profile file, and makes the
  92.  *        file more compact and easier to read. I have written a filter
  93.  *        'cvtprof' to convert from the old format to the new format, so
  94.  *        any existing profiles that people have can be easily converted.
  95.  * 
  96.  *        Profiles created with the -o option in interactive mode are written
  97.  *        in the new format.
  98.  * 
  99.  * Revision 1.8  90/11/16  16:33:45  16:33:45  dave (Dave Waller)
  100.  * Cleaned up rounding algorithm for depth profile; added round up to
  101.  * compartment bar graph display, so that it reads 100% properly;
  102.  * Fixed scale graphic at top of compartment bar graph (it was off
  103.  * by one character position to the right, contributing tothe bars not
  104.  * being at 100% visually when they actually were).
  105.  * 
  106.  * Revision 1.7  90/11/16  10:34:26  10:34:26  dave (Dave Waller)
  107.  * Added ability to specify different sample periods. Sample period ('s' in
  108.  * the source) is specified in the following manner (decreasing precedence):
  109.  * 
  110.  *            1) via command line option
  111.  *               Example: divecomp -mEdge -omyprof -s2.5
  112.  * 
  113.  *               Units are in minutes
  114.  * 
  115.  *            2) From within a profile file. The first line can optionally
  116.  *               specifiy a sample rate, if it has the form "s:<rate>".
  117.  *               Example: A profile that contains 2 minute samples would
  118.  *               look like
  119.  * 
  120.  *               s:2.0
  121.  *               60
  122.  *               60
  123.  *               60
  124.  *               60
  125.  *               .
  126.  *               .
  127.  *               .
  128.  * 
  129.  *               (the ellipses [...] are not part of the file). 'divecomp'
  130.  *               now automatically writes the sample period into every profile
  131.  *               it generates via the -o option.
  132.  * 
  133.  *             3) From a model file. Format is the same as a profile file.
  134.  *                The most basic command, "divecomp -mEdge", for instance,
  135.  *                will cause the computer to use the proper sampling rate
  136.  *                for the Edge. Also, any profiles generates with this model
  137.  *                will have the correct sampling value when used in different
  138.  *                models.
  139.  * 
  140.  *             4) Default: 1 minute.
  141.  * 
  142.  * Revision 1.6  90/11/14  17:36:52  17:36:52  dave (Dave Waller)
  143.  * Added compile time checks to modify curses behavior to accomodate
  144.  * incomplete curses library on Sun/OS. Controlling tissue is highlighted
  145.  * with a preceding '*' instead of reverse video, deiling message does not
  146.  * appear in reverse video when tissue M values are exceeded, and tissue
  147.  * bar graph is composed of '#' characters instead of reverse video
  148.  * spaces. Not quite as pretty as the full curses version, but that's what
  149.  * you get if you're running on a Sun.
  150.  * 
  151.  * Original graphical functionality is maintained for other platforms.
  152.  * 
  153.  * Revision 1.5  90/11/14  13:58:43  13:58:43  dave (Dave Waller)
  154.  * Added the following features:
  155.  * 
  156.  * - Decompression calculation and indication via ndl()
  157.  * - ingas/outgas indicator nxt to each compartment number
  158.  * - full "bottom timer" functions (i.e. dive #, bottom time, surface int)
  159.  * 
  160.  * Fixed the following bugs:
  161.  * 
  162.  * - would run forever if duration of 0 entered in interactive mode.
  163.  *   Program now rejects such an entry, and asks for depth and duration
  164.  *   again.
  165.  * 
  166.  * - Controlling tissue indication lagged one sample behind actual value;
  167.  *   Fixed.
  168.  * 
  169.  * - Dive profile graph now rounds up to the nearest depth increment.
  170.  * 
  171.  * Revision 1.4  90/11/13  16:38:49  16:38:49  dave (Dave Waller)
  172.  * Fixed up the comments.
  173.  * 
  174.  * Revision 1.3  90/11/13  15:11:48  15:11:48  dave (Dave Waller)
  175.  * Initial checkin; started using RCS to keep track of revisions. This
  176.  * revision also fixes a bug with the code that reads in half times from model
  177.  * files -- they were stored as type 'int', which created problems for
  178.  * fscanf if the half times were floats. Changed the array 'half[]' to
  179.  * type 'float'.
  180.  *  */
  181.  
  182. /* dive computer simulator, by Dave Waller (davew@hpdstma.hp.com)
  183.  
  184. ========================================================================
  185. Acknowlegments:
  186.  
  187. Much thanks to Eric Williams (sargon@portia.stanford.edu) for his
  188. outstanding contribution to the ndl() routine to incorporate ceilings
  189. and decompression times.  His support and timely reviews/bug reports
  190. have been extremely helpful.  Thanks, Eric!
  191. ========================================================================
  192.  
  193. This program simulates a multi-compartment model dive computer, based
  194. upon modified Haldanean tissue absoption models.  The number of
  195. compartments and their corresponding constants are not built into the
  196. program, but must be loaded at runtime, allowing the program to
  197. simulate most dive computers on the market today.
  198.  
  199. Subsequent derivation of equations as found in the comments below were
  200. done by myself, as I sit here at work, using only my memory;
  201. therefore, this initial pass at the simulation may contain some errors
  202. that I will have to fix later, after I have had time to review the
  203. pertinent material at home.  Surprisingly enough, it seems to work
  204. pretty accurately, based upon my experience with my ORCA Delpi.  Have
  205. fun!
  206.  
  207.        Tissue halftimes in general express the time it takes for a
  208.        tissue to either absorb or release nitrogen such that the
  209.        tissue is 50% equalized to the pressure differential between
  210.        the ambient nitrogen pressure and the tissue nitrogen
  211.        pressure. This uptake or outgassing obeys an asymptotic
  212.        exponential, namely if ambient pressure is PA, and initial tissue
  213.        pressure is PT0, then the tissue pressure as a function of
  214.        time is expressed as,
  215.  
  216.        PT = PT0 + (PA - PT0)(1 - e^(-kt))
  217.  
  218.        where k = -ln(0.5)/T     = 0.6931/T
  219.                     half           half
  220.  
  221.        This formula holds true for static conditions; i.e. PA does
  222.        not change. In reality, a diver will probably be continuously
  223.        changing depth, and therefore PA is changing also. A
  224.        reasonable approximation of continous PT values can be
  225.        calculated by sampling PA at discrete intervals, and applying
  226.        the above formula over the preceding sample interval, then
  227.        starting over with a new PA and PT0 for the next interval
  228.        based upon the calculation.
  229.  
  230.        In this simulated computer, we assume datapoints are recorded
  231.        one per minute. Since the halftimes are expressed in minutes,
  232.        the above formula can be reduce to an adjustment at each
  233.        sample that consists of the following:
  234.  
  235.        PTnew = PTlast + K(PA - PTlast)
  236.  
  237.        where K = (1 - e^(-k*1)) = 1-e^(-k) = 1 - (1/e^k) =
  238.  
  239.                    1 - (0.5)^(1/T    )
  240.                          half
  241.        
  242.        Note that the K value is dependant upon consistency between
  243.        the units of the compartment half time and the computer
  244.        sample period. Specifically, the above formula only applies
  245.        if the sample period is 1 unit of the half time unit (in this
  246.        case minutes). For an arbitrary sample period s, the formula
  247.        for the K value is
  248.  
  249.        K = 1 - (0.5)^(s/T    )
  250.                  half
  251.  
  252.        The dive computer simulator calculates the K values on
  253.        startup, with the sample period 's' being taken from the
  254.        following sources, in order of decreasing priority:
  255.  
  256.        1) a sample period specified on the command line with the
  257.           -s option
  258.  
  259.        2) A sample period specified in a profile file. This must
  260.           take precedence over a sample period in a model file,
  261.           otherwise the number of samples in the file will not
  262.           match the dive profile in real time. In order to compare
  263.           different models, the same sampling period must be used
  264.           for a single profile.
  265.  
  266.        3) A sample period specified in the model file
  267.  
  268.        4) The default sample period of 1 second
  269.  
  270.        */
  271.  
  272. char scale100[]=" 10| 20| 30| 40| 50| 60| 70| 80| 90|100|  % Allowed N2";
  273. char scale200[]=" 20| 40| 60| 80|100|120|140|160|180|200|  % Allowed N2";
  274.  
  275. #include <stdio.h>
  276. #include <curses.h>
  277.  
  278. #define NCOMP        24
  279. #define ONE_ATM_FSW    33.0
  280. #define PN2        0.79
  281. #define ONE_ATM_FSW_N2    PN2 * ONE_ATM_FSW
  282. #define LN_HALF        -0.6931
  283. #define INF         32000        /* for ndl() */
  284. #define DELTAM        1.0        /* temporary fix */
  285. #define SHOW_LOADING    0
  286. #define SHOW_FSWA    1
  287.  
  288. /* Some obscure variable names and their meanings:
  289.  
  290.     bt    Bottom Time
  291.  
  292.     si    Surface Interval
  293.  
  294.     cf    Compression factor (for compressing profile display)
  295.  
  296.     di    Depth Increment (for profile display)
  297.  
  298. */
  299.  
  300. int cf=1, duration=0, samples=0, n=0, controlling_tissue, lines=24,
  301.     ceiling, dive=1, mode=SHOW_LOADING, logmode=SHOW_LOADING,
  302.     decomode=0, width=80;
  303.  
  304. float bt=0.0, si=0.0, apn2=ONE_ATM_FSW_N2, half[NCOMP], ucf=1.0, di,
  305.     depth, Kvalues[NCOMP], Mvalues[NCOMP], c[NCOMP], s=0.0,
  306.     depth_record[8192];
  307.  
  308. double ceil(), floor(), pow(), exp(), log();
  309. float loading();
  310. char buf[128], unitstr[16], scale100FSWA[80], scale200FSWA[80];
  311.  
  312. char revision[]="$Header: divecomp.c,v 2.5 90/12/12 15:30:02 dave Exp $";
  313.  
  314. FILE *profile=NULL, *outfile=NULL, *logfile=NULL;
  315.  
  316. main(argc,argv)
  317. int argc;
  318. char *argv[];
  319. {    
  320.     FILE *model=NULL, *fopen();
  321.     int i, opt;
  322.     char *getenv(), *envstr, *cptr, *strcpy(), *strchr(), mname[32];
  323.     char *strcat();
  324.     float get_sample(), dummy;
  325.     double atof();
  326.     extern char *optarg;
  327.  
  328.     strcpy(unitstr,"ft");
  329.  
  330.     /* find out how many lines display has */
  331.  
  332.     if ((envstr=getenv("LINES"))!=0) lines=atoi(envstr);
  333.     if ((envstr=getenv("WIDTH"))!=0) width=atoi(envstr);
  334.  
  335.     /* process command line options */
  336.  
  337.     while ((opt=getopt(argc,argv,"p:m:o:s:u:l:P"))!=EOF) {
  338.         switch(opt) {
  339.             case 'p':
  340.                 if ((profile=fopen(optarg,"r"))==NULL) {
  341.                     perror(optarg);
  342.                     exit(1);
  343.                 }
  344.                 break;
  345.             case 'm':
  346.                 if ((model=fopen(optarg,"r"))==NULL) {
  347.                      perror(optarg);
  348.                      exit(1);
  349.                 }
  350.                 strcpy(mname,optarg);
  351.                 break;
  352.             case 'o':
  353.                 if ((outfile=fopen(optarg,"w"))==NULL) {
  354.                     perror(optarg);
  355.                     exit(1);
  356.                 }
  357.                 break;
  358.             case 's':
  359.                 s = atof(optarg);
  360.                 break;
  361.             case 'u':
  362.                 cptr = strchr(optarg,(int)':');
  363.                 *cptr = '\0';
  364.                 ucf = atof(cptr+1);
  365.                 strcpy(unitstr, optarg);
  366.                 break;
  367.             case 'l':
  368.                 if ((logfile=fopen(optarg,"w"))==NULL) {
  369.                     perror(optarg);
  370.                     exit(1);
  371.                 }
  372.                 break;
  373.             case 'P':
  374.                 mode=logmode=SHOW_FSWA;
  375.                 break;
  376.                 
  377.         }
  378.     }
  379.  
  380.     if (argc == 1) {
  381.         fprintf(stderr, "Valid options are:\n");
  382.         fprintf(stderr, " -m<fname>        Model file.  REQUIRED.\n");
  383.         fprintf(stderr, " -p<fname>        Input profile filename\n");
  384.         fprintf(stderr, " -o<fname>        Output filename for dive profile\n");
  385.         fprintf(stderr, " -l<fname>        Output logfile filename.  logfile is depth and N2 vs time.\n");
  386.         fprintf(stderr, " -s<float>        Sample period (in minutes)\n");
  387.         fprintf(stderr, " -u<name>:<float> Depth units.  (new-unit)*<float> = 1 ft\n");
  388.         fprintf(stderr, " -P               Display compartments as (units)SWA\n");
  389.         exit(0);
  390.     }
  391.  
  392.     if (model==NULL) {
  393.         fprintf(stderr, "You must specify a model with the \"-m\" option\n");
  394.         exit(0);
  395.     } 
  396.  
  397.     /* check profile for depth unit specification */
  398.  
  399.     if (profile)
  400.         if(fscanf(profile,"u:%[^:]:%f",unitstr,&ucf)==0)
  401.             rewind(profile);
  402.  
  403.     /* check model file for sampling period */
  404.  
  405.     if (fscanf(model,"s:%f",(s==0?&s:&dummy))==0) rewind(model);
  406.  
  407.     if (s == 0.0) s = 1.0; /* default */
  408.  
  409.     /* set up the pressure scales */
  410.  
  411.     for (i=1; i<11; i++) {
  412.         sprintf(buf,"%3d|",(int)(i*10*ucf));
  413.         strcat(scale100FSWA,buf);
  414.         sprintf(buf,"%3d|",(int)(i*20*ucf));
  415.         strcat(scale200FSWA,buf);
  416.     }
  417.  
  418.     sprintf(buf,"  %s sea water abs.", unitstr);
  419.     strcat(scale100FSWA,buf);
  420.     strcat(scale200FSWA,buf);
  421.  
  422.     /* read in the computer model */
  423.  
  424.     n = 0;
  425.     while(fscanf(model, "%f %f", &half[n], &Mvalues[n]) != EOF) {
  426.         Kvalues[n] = 1.0 - pow((double)0.5, (double)(s/half[n]));
  427.         n++;
  428.     }
  429.     fclose(model);
  430.  
  431.     if (profile && outfile) {
  432.         fprintf(stderr,"-o option can only be used in interactive mode.\n");
  433.         exit(1);
  434.     }
  435.  
  436.     /* put the depth unit into the outfile if there is one */
  437.  
  438.     if (outfile) fprintf(outfile,"u:%s:%.4f\n",unitstr,ucf);
  439.  
  440.     /* If logging, write the header into the logfile */
  441.  
  442.     if (logfile) {
  443.         fprintf(logfile,"Model: %s   Sample period: %f min\n\n",
  444.             mname, s);
  445.         if (logmode==SHOW_LOADING)
  446.             fprintf(logfile,"Depth\t%% Allowed Nitrogen\n(%s)\t",
  447.                 unitstr);
  448.         else fprintf(logfile,"Depth\tCompartment pressures (%s sea water abs.)\n(%s)\t", unitstr, unitstr);
  449.         for (i=0; i<n; i++) fprintf(logfile,"c%02d   ",i+1);
  450.         fprintf(logfile,"\n");
  451.     }
  452.  
  453.     /* determine the depth increment used in the profile display */
  454.  
  455.     for (di=ucf*5; ucf*200/di > lines-n-4; di+=ucf*5);
  456.  
  457.     /* initialize tissue N2 levels */
  458.  
  459.     for (i=0; i<NCOMP; i++) c[i] = ONE_ATM_FSW_N2;
  460.  
  461.     /* initialize the display */
  462.  
  463.     init_display();
  464.     update_display();
  465.  
  466.     for (;;) {
  467.  
  468.         get_sample();
  469.  
  470.         /* update time values */
  471.  
  472.         if (depth > 3) {
  473.  
  474.             /* all depths >3ft are counted as bottom time */
  475.  
  476.             if ((si > 0.0) && (si < 5.0)) {
  477.  
  478.                 /* We must account for brief periods of
  479.                 surfacing during the dive. The surface
  480.                 interval counter starts immediately upon
  481.                 ascending to depths shallower than 3
  482.                 feet, yet if the diver returns to depths
  483.                 deeper than 3 feet within 5 minutes, the
  484.                 time elapsed as surface interval must be
  485.                 added into the bottom time value, and
  486.                 the surface interval reset to 0 */
  487.  
  488.                 bt += si + s;
  489.                 si = 0.0;
  490.             }
  491.             else if (si == 0.0) bt+=s; /* no intermediate
  492.                            surface interval to
  493.                            account for */
  494.  
  495.             else if (si > 5.0) {       /* else the diver has
  496.                            been on the surface
  497.                            for more than 5
  498.                            minutes, so we treat
  499.                            this as a new dive
  500.                            */
  501.                 dive++;
  502.                 bt=s;
  503.                 si=0.0;
  504.             }
  505.         }
  506.         else si+=s; /* else the diver is still shallower than
  507.                 3 feet; update the surface interval counter
  508.                 */
  509.         samples++;
  510.  
  511.         /* update ambient N2 pressure */
  512.  
  513.         apn2 = PN2*(depth+ONE_ATM_FSW);
  514.  
  515.         /* update tissue (compartment) N2 levels */
  516.  
  517.         update_tissues();
  518.  
  519.         /* check if there has been a character input */
  520.  
  521.         switch (getch()) {
  522.             case 't':
  523.                 mode =
  524.                    (mode==SHOW_LOADING?SHOW_FSWA:SHOW_LOADING);
  525.                 break;
  526.  
  527.             case 'q':
  528.                 decomode = 0; /* reset vars that cause */
  529.                 duration = 0; /* continuous calculation */
  530.                 break;
  531.  
  532.             case KEY_DOWN:
  533.                 depth += di;
  534.                 break;
  535.  
  536.             case KEY_UP:
  537.                 depth -= di;
  538.                 break;
  539.         }
  540.  
  541.         /* repaint the display */
  542.  
  543.         update_display();
  544.         update_dive_profile();
  545.         refresh();
  546.  
  547.         /* if logging, write entry */
  548.  
  549.         if (logfile) {
  550.             fprintf(logfile,"%.1f\t", ucf*depth);
  551.             for (i=0; i<n; i++) {
  552.                 if (logmode == SHOW_LOADING)
  553.                     fprintf(logfile,"%-5d ",
  554.                         (int)(100*loading(i)));
  555.                 else fprintf(logfile,"%-5.1f ", c[i]);
  556.             }
  557.             fprintf(logfile,"\n");
  558.         }
  559.     }
  560. }
  561.  
  562. update_display()
  563. {
  564.     /* update the stats line first. Since the controlling tissue is
  565.     found in the ndl() routine which is called by the update_stats()
  566.     routine, it must be called before the tissue graph is repainted.
  567.     */
  568.  
  569.     update_stats();
  570.     update_ndl_limits();
  571.  
  572.     if (mode == SHOW_LOADING) display_loading();
  573.     else display_fswa();
  574. }
  575.  
  576. update_tissues()
  577. {
  578.     int i;
  579.  
  580.     for (i=0; i<n; i++) /* n == number of compartments */
  581.         c[i] = c[i] + (apn2 - c[i])*Kvalues[i];
  582. }
  583.  
  584. update_stats()
  585. {
  586.     int limit;
  587.     float time_to_surface();
  588.     char *showtime(), buf1[16], buf2[16], buf3[16];
  589.  
  590.     move(n+2,15);
  591.     clrtoeol();
  592.     limit = ndl(apn2);
  593.  
  594.     sprintf(buf, "Dive: %d  Depth: %.1f %s  BT: %s",
  595.         dive, ucf*depth, unitstr,  showtime(buf1,bt));
  596.     addstr(buf);
  597.     move(n+3,15);
  598.     clrtoeol();
  599.     if (limit >= 0) {
  600.         sprintf(buf,"NDL: %s  Sfc Int: %s",
  601.             showtime(buf2,(float)limit), showtime(buf3,si));
  602.         addstr(buf);
  603.     }
  604.     else {
  605.         attron(A_REVERSE);
  606.         sprintf(buf,"Ceiling: %.1f %s", ucf*ceiling,
  607.             unitstr);
  608.         addstr(buf);
  609.         sprintf(buf," Time to surface: %s",
  610.             showtime(buf1,time_to_surface()));
  611.         addstr(buf);
  612.         attroff(A_REVERSE);
  613.     }
  614. }
  615.  
  616. update_dive_profile()
  617. {
  618.     int i, j, x, y;
  619.  
  620.     /* update the dive profile display */
  621.  
  622.     if (samples/cf > 60) {
  623.  
  624.         /* erase and compress display */
  625.  
  626.         for (i=n+5; i<lines; i++) {
  627.             move(i,10);
  628.             clrtoeol();
  629.         }
  630.  
  631.         cf *=2;
  632.         for (i=1; i<samples; i+=cf) {
  633.             x = 10+i/cf;
  634.             y = n+5+(int)ceil((double)ucf*depth_record[i]/di);
  635.             move((y>=lines?lines-1:y),x);
  636.             addch((y>=lines?'v':'*'));
  637.         }
  638.     }
  639.  
  640.     x = 10+samples/cf;
  641.     y = n+5+ceil((double)ucf*depth/di);
  642.     move((y>=lines?lines-1:y),x);
  643.     addch((y>=lines?'v':'*'));
  644.  
  645.     depth_record[samples] = depth;
  646. }
  647.  
  648. update_ndl_limits()
  649. {
  650.     int i;
  651.     int t;
  652.     float d;
  653.     float p;
  654.     char buf[10];
  655.  
  656.     for (i=n+5; i<lines; i++) {
  657.         d = (i-n-5)*di/ucf;        /* depth in feet */
  658.         p = 0.79 * (d + 33.0);        /* pn2 at depth  */
  659.         t = ndl(p);            /* compute limits */
  660.         if (t<0) t=0;
  661.  
  662.         move(i,0);
  663.         if (t >= INF) {
  664.             addstr("INF");
  665.         } else if (t >= 1000) {
  666.             addstr("***");
  667.         } else if (t>0) {
  668.             sprintf(buf,"%3d",t);
  669.             addstr(buf);
  670.         } else addstr("   ");
  671.     }
  672.     t = ndl(apn2);    /* leave all the globals at the right values */
  673. }
  674.  
  675. display_loading()
  676. {
  677.     int i, j, k, temp, over;
  678.     /* if overpressure on any tissue, change the scale to 0 - 200% */
  679.  
  680.     for (over=i=0; i<n; i++) if (c[i] > Mvalues[i]) over=1;
  681.     move(0,10);
  682.     clrtoeol();
  683.     if (over) addstr(scale200);
  684.     else addstr(scale100);
  685.  
  686.     /* display tissue saturation */
  687.  
  688.     for (i=0; i<n; i++) { /* n == number of compartments */
  689.  
  690.         /* display the compartment number, highlighting if it
  691.         is the controlling tissue. */
  692.  
  693. #ifdef SIMPLE
  694.         move(1+i,6);
  695.         if (i==controlling_tissue) sprintf(buf,"*%d",i+1);
  696.         else sprintf(buf," %d",i+1);
  697. #else
  698.         move(1+i,7);
  699.         if (i==controlling_tissue) attron(A_REVERSE);
  700.         sprintf(buf,"%d",i+1);
  701. #endif
  702.         addstr(buf);
  703. #ifndef SIMPLE
  704.         attroff(A_REVERSE);
  705. #endif
  706.  
  707.         /* display the ingas/outgas indicator for the
  708.         compartment. Ingas/outgas indicator is adjusted to 3 decimal
  709.         place accuracy. */
  710.  
  711.         move(1+i,9);
  712.         temp = (apn2 - c[i]) * 1000.0;
  713.         if (temp > 0) addch('>');
  714.         if (temp == 0) addch(' ');
  715.         if (temp < 0) addch('<');
  716.  
  717.         /* clear the bar before displaying the new value */
  718.  
  719.         move(1+i,10);
  720.         clrtoeol();
  721.  
  722.         k = ceil(floor((double)(1000.0 * loading(i))) * (over==0?0.04:0.02));
  723.  
  724. #ifndef SIMPLE
  725.         attron(A_REVERSE | A_UNDERLINE);
  726. #endif
  727.         sprintf(buf," %4.2f",ucf*c[i]);
  728.         for (j=0; j<k; j++) {
  729.             if (j==(width-11-strlen(buf))) {
  730.                 addch('>');
  731.                 break;
  732.             }    
  733. #ifdef SIMPLE
  734.             if (j<(over==0?40:20)) addch('#');
  735. #else
  736.             if (j<(over==0?40:20)) addch(' ');
  737. #endif
  738.             else addch('*');
  739.         }
  740. #ifndef SIMPLE
  741.         attroff(A_REVERSE | A_UNDERLINE);
  742. #endif
  743.         addstr(buf);
  744.     }
  745. }
  746.  
  747.  
  748. display_fswa()
  749. {
  750.     int i, j, k, temp, over, Mmarker;
  751.     /* if overpressure on any tissue, change the scale to 0 - 200% */
  752.  
  753.     for (over=i=0; i<n; i++) if (c[i] > 100.0) over=1;
  754.     move(0,10);
  755.     clrtoeol();
  756.     if (over) addstr(scale200FSWA);
  757.     else addstr(scale100FSWA);
  758.  
  759.     /* display compartment N2 levels in FSWA */
  760.  
  761.     for (i=0; i<n; i++) { /* n == number of compartments */
  762.  
  763.         /* display the compartment number, highlighting if it
  764.         is the controlling tissue. */
  765.  
  766. #ifdef SIMPLE
  767.         move(1+i,6);
  768.         if (i==controlling_tissue) sprintf(buf,"*%d",i+1);
  769.         else sprintf(buf," %d",i+1);
  770. #else
  771.         move(1+i,7);
  772.         if (i==controlling_tissue) attron(A_REVERSE);
  773.         sprintf(buf,"%d",i+1);
  774. #endif
  775.         addstr(buf);
  776. #ifndef SIMPLE
  777.         attroff(A_REVERSE);
  778. #endif
  779.  
  780.         /* display the ingas/outgas indicator for the
  781.         compartment. Ingas/outgas indicator is adjusted to 3 decimal
  782.         place accuracy. */
  783.  
  784.         move(1+i,9);
  785.         temp = (apn2 - c[i]) * 1000.0;
  786.         if (temp > 0) addch('>');
  787.         if (temp == 0) addch(' ');
  788.         if (temp < 0) addch('<');
  789.  
  790.         /* clear the bar before displaying the new value */
  791.  
  792.         move(1+i,10);
  793.         clrtoeol();
  794.  
  795.         Mmarker = ceil((double)(Mvalues[i]*(over==0?0.4:0.2)));
  796.         k = ceil((double)(c[i]*(over==0?0.4:0.2)));
  797.             
  798.  
  799. #ifndef SIMPLE
  800.         attron(A_REVERSE | A_UNDERLINE);
  801. #endif
  802.         for (j=0; j<k; j++) {
  803.             if (j==(width-11)) {
  804.                 addch('>');
  805.                 break;
  806.             }    
  807. #ifdef SIMPLE
  808.             if (j<Mmarker) addch('#');
  809. #else
  810.             if (j<Mmarker) addch(' ');
  811. #endif
  812.             else addch('*');
  813.         }
  814. #ifndef SIMPLE
  815.         attroff(A_REVERSE | A_UNDERLINE);
  816. #endif
  817.     }
  818. }
  819.  
  820. init_display()
  821. {
  822.     int i;
  823.     initscr();
  824.  
  825.     /* set nodelay mode for input processing */
  826.  
  827.     nodelay(stdscr,TRUE);
  828.     keypad(stdscr,TRUE);
  829.  
  830.     /* print the saturation %age scale */
  831.  
  832.     move(0,10);
  833.     if (mode==SHOW_LOADING) addstr(scale100);
  834.     else addstr(scale100FSWA);
  835.  
  836.     /* label the Dive profile graph */
  837.  
  838.     move(n+3,0);
  839.     addstr("NDL  Depth");
  840.     move(n+4,0);
  841.     addstr("----------");
  842.  
  843.     /* create the dive profile display */
  844.  
  845.     for (i=n+5; i<lines; i++) {
  846.         move(i,4);
  847.         sprintf(buf,"%5.1f|",(i-n-5)*di);
  848.         addstr(buf);
  849.     }
  850. }
  851.  
  852. float get_sample()
  853. {
  854.     float target_depth;
  855.     double atof();
  856.     while (duration == 0) {
  857.         if (!profile) {
  858.  
  859.             /* check for autodecompression mode */
  860.  
  861.             if (decomode) {
  862.                 if (ceiling == 0) decomode = 0;
  863.                 else {
  864.                     target_depth =
  865.                         ceil((double)ceiling/5.0)*5.0;
  866.                     if (target_depth - depth > 30)
  867.                         depth -= 30;
  868.                     else depth = target_depth;
  869.                     duration = 1;
  870.                 }    
  871.                 continue;
  872.             }
  873.  
  874.             /* read depth and duration from user */
  875.  
  876.             nodelay(stdscr,FALSE);
  877.             move(n+4,15);
  878.             clrtoeol();
  879.             addstr("Enter Depth ('q' to quit): ");
  880.             refresh();
  881.             getstr(buf);
  882.             if (buf[0]=='q') {
  883.                 nodelay(stdscr,FALSE);
  884.                 endwin();
  885.                 exit(0);
  886.             }
  887.  
  888.             /* 'd' indicates autodecompression */
  889.  
  890.             if (buf[0]=='d') {
  891.                 decomode = 1;
  892.                 move(n+4,15);
  893.                 clrtoeol();
  894.                 attron(A_REVERSE);
  895.                 addstr("*** autodecompressing ***");
  896.                 attroff(A_REVERSE);
  897.                 continue;
  898.             }
  899.  
  900.             depth = atof(buf);
  901.             move(n+4,15);
  902.             clrtoeol();
  903.             addstr("Enter Duration (min): ");
  904.             refresh();
  905.             getstr(buf);
  906.             duration = atof(buf) / s;
  907.             move(n+4,15);
  908.             clrtoeol();
  909.             if (outfile && duration)
  910.                 fprintf(outfile,"%.1f %d\n",depth,
  911.                     (int)(duration*s));
  912.         }
  913.  
  914.         /* else get depth and duration from profile */
  915.  
  916.         else if (fscanf(profile,"%f %d",&depth, &duration)==EOF) {
  917.             endwin();
  918.             exit(0);
  919.         }
  920.         else duration /= s;
  921.         depth /= ucf;
  922.     }
  923.     duration--;
  924.     nodelay(stdscr,TRUE);
  925.     return(depth);
  926. }
  927.  
  928. float time_to_surface()
  929. {
  930.     int i, j;
  931.     float time=0, ctemp[NCOMP], stop, next, pn2_stop, pn2_next,
  932.           target, temp, dt;
  933.  
  934.     /* make a copy of the compartment N2 partial pressure values */
  935.  
  936.     for (i=0; i<n; i++) ctemp[i] = c[i];
  937.     
  938.     /* set up first stop data */
  939.  
  940.     stop = ceil((double)(ceiling/5.0))*5.0;
  941.     pn2_stop = PN2*(stop+ONE_ATM_FSW);
  942.  
  943.     /* for each stop, calculate the time to get to it from the
  944.     previous stop */
  945.  
  946.     for (next=stop-5.0; next>=0.0; next-=5.0) {
  947.         pn2_next = PN2*(next+ONE_ATM_FSW);
  948.         for (i=0, temp=0; i<n; i++) {
  949.             /* determine the time to decompress compartment
  950.             to target N2 pressure. Target pressure is the
  951.             value at which the compartment yields a ceiling
  952.             equal to the next stop. Time is then calculated
  953.             using the exponential decay formula with the
  954.             current stop as the ambient pressure */
  955.  
  956.             target = next * DELTAM + Mvalues[i];
  957.             if ((target < ctemp[i])&&(pn2_stop < ctemp[i])) {
  958.                 dt = half[i] / LN_HALF *
  959.                     log((double)(1.0 -
  960.                     (target - ctemp[i]) /
  961.                     (pn2_stop - ctemp[i])));
  962.                 if (dt > temp) temp = dt;
  963.             }
  964.  
  965.         }
  966.  
  967.         /* convert dt to a number of samples */
  968.  
  969.         dt = ceil((double)(temp/s));
  970.  
  971.         /* adjust compartment values for next calculation */
  972.  
  973.         for (i=0; i<n; i++)
  974.             for (j=0; j<dt; j++) ctemp[i] = ctemp[i] +
  975.                 (pn2_stop - ctemp[i]) * Kvalues[i];
  976.         time += dt*s;
  977.         stop = next;
  978.         pn2_stop = PN2*(stop+ONE_ATM_FSW);
  979.     }
  980.     return(time);
  981. }
  982.  
  983. ndl(pressure)
  984. float pressure;
  985. {
  986.     int i, limit, dlimit;
  987.     int dcontrol;
  988.     int safe, outgas, decom;
  989.     int t;
  990.     double time_fn();
  991.  
  992.     limit = INF;
  993.     dlimit = 0;
  994.     dcontrol = -1;
  995.     ceiling = 0;
  996.     for (i=0; i<n; i++) {
  997.         safe = (pressure < Mvalues[i]);
  998.         outgas = (pressure < c[i]);
  999.         decom = (Mvalues[i] < c[i]);
  1000.         if (decom) {
  1001.             /* compute ceiling */
  1002.             t = ceil((c[i] - Mvalues[i])/DELTAM);
  1003.             if (t > ceiling) ceiling = t;
  1004.  
  1005.             if ((outgas) && (safe)) {
  1006.                 /* how long to decompress */
  1007.                 t = - ceil(time_fn(pressure,i));
  1008.                  if (t < dlimit) { 
  1009.                     dlimit = t;
  1010.                     dcontrol = i;
  1011.                 }
  1012.             } else {
  1013.                 /* compartment will not decompress */
  1014.                 dlimit = -INF;
  1015.                 dcontrol = i;
  1016.             }
  1017.             } else if ((!safe) && (!outgas)) {
  1018.              /* compute no-decompress time */
  1019.  
  1020.             if (pressure == Mvalues[i]) {
  1021.                 t = INF;
  1022.             } else if (pressure == c[i]) {
  1023.                 t = INF;
  1024.             } else if (c[i] == Mvalues[i]) {
  1025.                 t = 0;
  1026.             } else {
  1027.                 t = floor(time_fn(pressure,i));
  1028.             }
  1029.  
  1030.             if (t < limit) {
  1031.                 limit = t;
  1032.                 controlling_tissue = i;
  1033.             }
  1034.         }
  1035.     }
  1036.     
  1037.     /* if any tissues are over their m-value than they will control */
  1038.     if (dcontrol >=0) {
  1039.         controlling_tissue = dcontrol;
  1040.             limit = dlimit;
  1041.     }
  1042.           
  1043.     return(limit);
  1044. }
  1045.  
  1046. char *showtime(tmpbuf,min)
  1047. char *tmpbuf;
  1048. float min;
  1049. {
  1050.     char *strcpy();
  1051.     if ((int)min >= INF) return(strcpy(tmpbuf,"INF"));
  1052.     if ((int)min > 60) sprintf(tmpbuf,"%d h %d m", (int)min/60, 
  1053.         (int)min%60);
  1054.     else sprintf(tmpbuf,"%d min", (int)min);
  1055.     return(tmpbuf);
  1056. }
  1057.  
  1058. double time_fn(pressure, i)
  1059. float pressure;
  1060. int i;
  1061. {
  1062.     double temp;
  1063.  
  1064.     temp = (pressure-Mvalues[i])/(pressure-c[i]);
  1065.     return(half[i]/LN_HALF * log(temp));
  1066. }
  1067.  
  1068. float loading(index)
  1069. int index;
  1070. {
  1071.     return((c[index]-ONE_ATM_FSW_N2)/(Mvalues[index]-ONE_ATM_FSW_N2));
  1072. }
  1073.